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1.  Introduction 

We  consider  the  problem  of  solving  an  elliptic  partial  differential  equation  on  a  domain  Q  that 
is  broken  up  into  subregions  ft,.  By  using  domain  decomposition  or  substructuring  techniques,  the 
problem  reduces  to  separately  solving  approximate  problems  in  the  subdomains  and  updating  the 
solution  at  the  interfaces  between  two  or  more  subregions.  There  are  several  reasons  why  these 
techniques  might  be  attractive: 


1.  Special  solution  techniques,  like  fast  direct  solvers,  might  exist  to  solve  the  problems  on  the 
subdomains  that  cannot  be  applied  efficiently  to  the  entire  domain.  This  is  often  the  case,  for 
example,  when  the  subdomains  fi,  have  very  regular  geometry,  but  Q  does  not. 

2.  The  equations  in  the  different  subdomains  might  have  different  parameters  or  even  be  of  dif¬ 
ferent  nature,  in  which  case  the  idea  of  substructuring  comes  very  naturally. 

3.  The  idea  is  attractive  for  parallel  processing,  since  the  problem  can  be  decoupled  into  in¬ 
dependent  subproblems  and  the  communication  needed  is  limited  to  the  boundaries  of  the 
subdomains. 


For  the  class  of  domain  decomposition  methods  considered  in  this  paper,  the  basic  idea  consist 
of  the  following:  the  domain  is  discretized  and  partitioned  into  several  subregions,  then,  by  ap¬ 
plying  block  elimination  to  the  discretized  equations,  a  system  is  derived  for  the  unknowns  on  the 
interfaces  between  subregions.  This  system  is  sometimes  called  the  capacitance  system.  Forming 
the  right  hand  side  for  the  interface  system  and  solving  it  requires  the  solution  of  independent 
elliptic  problems  on  the  subdomains.  For  certain  constant  coefficient  problems,  fast  direct  methods 
can  be  applied  to  the  solution  of  the  interface  system.  Such  is  not  the  case,  however,  for  more 
general  operators  on  irregular  domains.  For  efficiency  reasons  the  system  must  then  be  solved  by 
iterative  methods,  such  as  the  preconditioned  conjugate  gradient  method  (PCG).  Once  the  solution 
is  known  on  the  interfaces,  one  more  elliptic  problem  must  be  solved  on  each  subdomain  with  the 
computed  values  as  boundary  conditions.  The  method  is  particularly  suited  to  problems  for  which 
the  subproblems  can  be  solved  efficiently,  for  example,  when  the  operator  has  separable  coefficients 
and  the  domain  is  a  union  of  rectangles.  On  the  other  hand,  when  the  subdomain  problems  cannot 
be  solved  efficiently  but  they  can  be  approximated  by  simpler  operators,  it  is  possible  to  derive 
block  preconditioners  for  the  original  system  based  on  preconditioners  for  the  capacitance  matrix. 

In  section  2.  we  illustrate  the  method  for  the  case  of  a  domain  that  is  the  union  of  two 
rectangles.  In  section  3.  we  consider  the  Poisson  and  Helmholtz  equations  on  a  rectangular  domain 
divided  into  parallel  strips  and  derive  the  capacitance  system  for  the  interface  variables.  For  these 
simple  and  regular  cases,  the  capacitance  system  can  be  solved  by  fast  direct  methods.  Such  is  not 
the  case  for  irregular  domains.  In  section  4.  we  summarize  the  various  preconditioners  that  were 
proposed  in  the  literature  for  use  with  the  CG  method.  For  the  case  of  variable  coefficient  problems, 
when  fast  direct  methods  are  not  applicable  to  the  solution  of  the  problems  on  the  subdomains, 
the  system  for  the  whole  domain  must  be  solved  by  iterative  methods.  Using  the  results  of  the 


previous  sections,  preconditioners  for  the  large  system  can  be  derived  from  preconditioners  for  the 
capacitance  matrix.  We  discuss  this  case  in  sectiou  5.  Finally,  in  section  6.  we  propose  a  new  family 
of  row-sum  preserving  banded  preconditioners  for  the  capacitance  matrix.  These  preconditioners 
have  the  advantage  that  they  can  be  applied  to  a  more  general  class  of  problems,  since  as  opposed 
to  most  of  the  other  preconditioners.  they  do  not  depend  on  special  properties  of  the  differential 
operator. 


2.  Domain  Decomposition 

We  will  first  consider  the  problem: 


La  =  f 


on  n 


(2.1)  'S 


Figure  1:  The  domain  fl  and  its  partition. 


with  boundary  conditions 

u  =  14  on  OCl 

where  I  is  a  linear  elliptic  operator  and  the  domain  Cl  is  as  illustrated  in  Fig.  1  .  We  will  call  the 
interface  between  Hi  and  CI2,  I\ 

If  we  order  the  unknowns  for  the  internal  points  of  the  subdomains  first  and  then  those  in  the 
interface  I\  then  the  discrete  solution  vector  u  ~  (14.  1/2.  U3)  satisfies  the  linear  system 


An  —  b  . 

that  can  be  expressed  in  block  form  as: 

( -An  -A13  \  f  1(1  \  /  61  \ 

-A22  -A23  1  "2  I  =  ^2  I 

V  -An  -A«  -A33  j  \  »3  )  V  ) 


(2.2) 


(2.3) 


The  system  (2.3)  can  be  solved  by  Block-Gaus«ian  Elimination  as  follows: 
Step  1:  Compute 

C  =  .A;w  -  -A^.AjY  .4|.»  -  -An-Aj_,1 


and  solve 


"1  —  -An% 

»■•>  =  A 22  b-2 

Cii.i  =  h  t  -  a[<  ir,  -  Al;' '>  ■.< 


(2.4) 

(2.5) 

(2.6) 


(2.7) 


Step  2:  Compute 


Note  that,  except  for  (2.7),  the  algorithm  only  requires  the  solution  of  problems  with  .4n  and 
.422,  which  corresponds  to  solving  independent  problems  on  the  subdomains.  This  technique  of 
reducing  the  problem  on  ft  to  the  solution  of  decoupled  problems  on  the  subdomains  and  a  smaller 
system  for  the  interface  is  usually  called  domain  decomposition  or  substructurini).  The  matrix  C 
(2.4)  is  the  Schur  complement  of  .433  in  .4  and  it  is  sometimes  called  the  capacitance  matrix  in  this 
context. 

3.  Poisson  and  Helmholtz  Equations  on  a  rectangle 

We  now  consider  the  case  where  L  is  the  Lapiaciau  operator  and  ft  is  a  rectangle  divided  into 
two  or  more  strips  like  is  shown  in  Fig.  2.  For  this  case,  the  exact  eigenvectors  and  eigenvalues  of 
C  are  known  [2,  6,  7]. 
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Figure  2:  Rectangular  domain  divided  into  strips. 

For  the  case  of  two  strips.  C  has  the  following  eigenvalue  decomposition: 
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t -T 


A  J 


where  IT'  is  the  matrix  whose  columns  are 


V  «  4-1 


(sin  jr/i.sin 2/-/i.  ■  .sin  njxh) 


and 


({- +  i  1  +  -  "'-+l  \  /  rr- 

XJ  =  ~  T - 7777T  +  7 - 77TTT  \  ni  *  -f- 


(3.1) 


(3.2) 


(3.3) 


for  j  =  1 , . . . ,  n  ,  where 


h  is  the  grid  size,  and  mi  and  m2  are  the  number  of  rows  of  grid  points  in  the  y-direction  in  f and 
CI2  respectively.  By  using  the  decomposition  (3.1),  the  capacitance  system  (2.7)  can  be  solved  by 
fast  Fourier  transforms.  Once  the  solution  m3  on  the  interface  is  computed,  we  can  compute  «i  and 
t/2  by  (2.8)  and  (2.9),  which  correspond  to  solving  two  independent  problems  on  the  subdomains 
with  boundary  condition  U3  on  T. 

In  the  multistrip  case,  the  matrix  C  has  the  block-tridiagonal  structure: 


/ Ci  B2  \ 

Bo  C2 

••  Bk 

V  Bk  Ck  J 


(3.6) 


all  blocks  C,  and  B,  have  the  same  matrix  of  eigenvectors  IV ,  i.e.  for  i  =  1  ,...,£,  we  have 

UrC.U'  =  A,  =  diag(A(1 . Ain)  (3.7) 

and  for  i  =  2, ....  k .  we  have 

=  D,  =  diag^S,, . 6in)  (3.8) 

where 


The  capacitance  system  can  then  be  solved  by  fast  Fourier  tran>forms  and  the  solution  of  n  decou¬ 
pled  tridiagonal  systems  of  dimension  /.-.  where  k  1  is  the  number  of  subdomains  (7j. 

Although  it  first  appears  that  the  algorithm  requires  the  solution  of  two  problems  on  each 
subdomain,  one  for  computing  the  right  hand  side  and  one  for  computing  the  solution  on  each 
subdomain,  some  computations  can  be  saved.  We  refer  the  interested  reader  to  [S],  where  a  detailed 
operation  count  is  derived  for  the  sequential  and  parallel  implementations. 

Formulas  (3.6)  to  (3.10)  can  be  generalized  to  two  particular  operators  other  than  the  Lapla- 
cian:  the  linear  elliptic  operator 

"  ri  +  •  (3.1 1) 

where  the  coefficient  j  take*  constant  values  3,  on  each  subdomain  F2,  and  the  Helmholtz  operator 

+  .  (3.12) 


4 


The  capacitance  matrix  for  the  operator  (3.11)  has  the  same  form  as  (3.6).  The  eigenvalues  of  C, 
and  B,  are  given  by 


and 


where 


1  +  »i+i.> 


1  -  u+ij 


+  3,+iVj 


(3.13) 


(3.14) 


The  operator  (3.11)  can  be  used  as  a  preconditioner  for  more  general  variable  coefficient  problems. 

The  Helmholtz  operator  (3.12)  has  important  applications  in  the  solution  of  time  dependent 
problems.  The  capacitance  matrix  for  this  operator  also  has  the  form  (3.6).  where  the  eigenvalues 
of  C,  and  B,  are  _ 


A 


>} 


vi-'/)m,+1  +  i-7;n‘+i+1J 


(3.15) 


and 


where 


(3.16) 

(3.17) 


and 


(3.18) 


4.  Poisson  Equation  on  Irregular  Domains 

In  general,  when  fi  has  irregular  shape  like  in  Fig.  1.  the  eigenvalues  and  eigenvectors  of  the 
capacitance  matrix  are  not  known.  The  computation  of  the  capacitance  matrix  is  expensive,  since 
it  requires  the  solution  of  m  +  1  systems  with  .4h  and  -djo-  and  it  is  also  expensive  to  invert  for  m 
large,  because  it  is  dense. 

Instead  of  solving  the  system  (2.7)  directly,  preconditioned  conjugate  gradient  methods  (PCG) 
can  be  applied,  where  only  matrix  vector  products  Cy  for  arbitrary  y  6  ltm  are  required.  This 
product  can  be  computed  by  solving  one  Poisson  problem  on  each  subdomain  with  boundary 
condition  on  T  given  by  y. 

Since  eacli  iteration  involves  the  solution  of  problems  on  the  subdomains,  keeping  t lie  number 
of  iterations  small  is  very  important  for  the  efficiency  of  the  method.  This  can  be  achieved  by 
choosing  a  good  preconditioner  for  C .  several  of  which  are  given  in  t In*  literature  (3.  10.  13.  5],  We 
summarize  these  here: 

1.  In  ^10'.  Dryja  proposes 

Md  =  vTf  . 


as  a  preconditioner  for  (2.4).  where  K  is  the  one-dimensional  Laplacian.  He  proved  that  the 
condition  number  of  the  preconditioned  system.  K(Mq1C)  is  bounded  independently  of  the 
mesh  size  h.  Since  Md  has  the  following  factorization 


Md  =  -.X»)\VT  . 

(4.1) 

where  the  columns  of  IV  are  given  by  (3.2)  and 

Xf  =  -2^/cTj 

(4.2) 

with  ct j  given  by  (3.4).  Md  can  be  inverted  by  FFT's. 

2.  Golub  and  Mayers  [13]  propose  a  preconditioner  given  by 

Mg  =  \f A"2  +  4/v  , 

which  has  the  following  decomposition: 

Mg  =  Wdiag(  .A£)irT  . 

(4.3) 

where 

III 

1 

10 

+ 

(4.4) 

Empirical  results  in  [13]  show  that  Mrj  performs  better  than  Md- 
3.  Another  interesting  preconditioner  was  given  by  Bjorstad  and  Widlund  [3]  and  has  the  following 
form: 

-'fa  =  -4:w  -  2.4f3.4j"11.4i3  . 

It  is  easy  to  show  that  the  eigenvalue  decomposition  of  Mb  is 

MB  =  Wdiag(  Af.A.f.  .  (4.5) 


where 


for  j  —  1 . n.  When  fii  and  n?  are  identical.  Mb  is  an  exact  preconditioner.  To  implement 

the  method.  Bjorstad  and  Widlund  solve  a  mixed  Neumann-Dirichlet  problem  in  one  of  the 
subdomains  and  a  Dirichlet  problem  in  the  other  one.  Their  method  has  the  advantage  that 
it  can  be  applied  to  more  general  operators  and  domain  shapes,  but  in  the  particular  case  of 
the  Laplacian  operator  on  a  union  of  rectangles,  it  is  less  efficient  than  applying  a  single  FFT 
computation  on  the  interface  grid  points,  as  the  factorization  (4.5)  sugests. 

4.  Although  Mp-  -Vc;  and  Mb  were  derived  independently  of  the  factorization  (3.1).  they  can  be 
viewed  as  progressively  better  approximations  to  the  capacitance  matrix  C  1 .  The  factorization 
(3.1)  is  exact  for  the  case  of  a  rectangular  fl.  while  Md  and  Mrj  are  not.  It  can  be  easily  observed 
that  (4.2)  is  a  first  order  approximation  to  (3.3).  while  (4.4)  is  a  second  order  approximation. 
On  the  other  hand.  Mg  i*  exact  only  for  the  case  of  a  rectangular  domain  divided  into  two 


1  Anderson  [l1  gives  an  interpretar ion  of  the  various  discrete  precon* lit  inner**  a*  approximations  to  a  continuous  operator  on 
the  interface 


identical  subregions.  All  this  sugests  that  (3.1)  might  be  a  better  preconditoner  for  the  case  of 
an  irregular  domain  [6]..  We  will  call  this  preconditioner  Me¬ 
in  Fig.  3  we  compare  the  preconditioners  Mp,  Mg  and  Me  for  the  Poisson  equation  on  a 
T-shaped  region  Cl  as  given  in  Fig.  1,  where  we  vary  the  aspect  ratio  of  the  subdomain  fii.  We 
consider  a  uniform  grid  on  Cl  with  15  grid  points  on  the  interface  T.  By  varying  mi,  the  number 
of  interior  grid  points  in  the  y  direction  on  the  subdomain  fti,  we  computed  the  condition  number 
of  the  preconditioned  capacitance  system  for  different  aspect  ratios  As  we  can  see,  Me 

performs  very  well,  even  when  fti  becomes  very  narrow,  while  the  others  deteriorate.  All  Me,  Me 
and  M\\r  are  indistinguishable  for  aspect  ratios  larger  than  one  and  they  are  all  better  than  Mp.  In 
[6],  Chan  analyzes  and  compares  these  preconditioners  on  rectangular  regions.  By  his  analysis,  we 
can  see  that  Me  is  always  better  than  Mq  on  a  rectangle  and  both  preconditoners  perform  poorly 
when  the  aspect  ratio  for  the  dimension  of  the  rectangles  is  small.  See  [14]  for  a  careful  numerical 
comparison  of  these  and  other  preconditioners  for  constant  and  variable  coefficients  operators. 


o«p«oi  roilo 


Figure  3:  T-shaped  region.  Condition  number  of  the  pre¬ 
conditioned  capacitance  matrix  with  Chans  (C),  Dryja's(D), 
Bjorstad  and  Widluud’s  (W)  and  Golub  and  Mayers’  (G)  pre¬ 
conditoners. 


5.  Variable  Coefficient  Problems 

In  the  case  of  non-constant  coefficient  problems,  there  usually  are  no  fast  solvers  available  for 
.4,,  an  .  A22  and  therefore  the  solutions  to  systems  with  these  matrices  are  to  be  avoided.  In  that 
case,  a  Krylov  subspace  method  can  be  applied  to  solve  the  system  (2.3)  on  the  whole  domain 
instead  of  just  the  capacitance  system  on  the  interface.  Therefore,  we  must  now  be  concerned  with 
the  problem  of  finding  preconditioners  for  (2.3)  that  make  use  of  the  regularity  of  the  subdomains. 
We  will  show  that  preconditioners  for  (2.3)  can  be  derived  from  preconditioners  for  the  capacitance 
matrix.  .Assume  that  the  variable  coefficient  operator  can  be  approximated  by  operators  with 
constant  coefficients  on  each  subdomain.  In  particular,  let  B\\  and  S2 2  be  approximations  to  .4,, 
and  .422,  corresponding  to  the  discretization  of  linear  elliptic  operators  with  constant  coefficients 
on  Q,  and  fii  respectively.  Based  on  the  following  decomposition  of  the  matrix  .4  in  (2.3): 


.4  = 


-4,i 

\  1 

.4 

.422 

/  .4-, 

-4.31 

-4.32  Cl  \ 

/ 

(5.1) 


where  C  is  the  Schur  complement  (2.4),  we  can  derive  a  preconditoner  for  .4  given  by: 


M  = 


Bh  -4 13 
^22*  *^23 


(5.2) 


where  M  is  a  good  preconditioner  for  the  matrix  C.  We  can  see  that  M  is  easily  invertible  by 
block-elimination,  since  fast  solvers  can  be  applied  to  solve  systems  with  B ,,  and  B22- 

Preconditioners  of  the  form  (5.2)  were  first  used  by  Bramble,  Pasciak  and  Schatz's  [4,  5).  In 
[4],  Drvja's  preconditioner  is  used  as  the  matrix  M  in  (5.2).  The  second  preconditioner  in  [5] 
corresponds  to  chosing  the  matrix  M  given  by  Bjbrst.ad  and  Widluud  [3).  As  a  generalization  of 
their  idea,  any  of  the  preconditioners  given  for  the  coustant  coefficients  case  can  be  applied  here  as 
M.  In  fact,  a  theorem  by  Eisenstat  in  [14]  shows  that,  when  B„  =  .4„,  the  PCG  algorithm  applied 
to  (2.7)  with  preconditioner  M  aud  initial  guess  k®  is  equivalent  to  the  PCG  algorithm  applied  to 
(2.3)  with  preconditioner  given  by  (5.2)  and  initial  guess  ( .4 j-,1  (61  -  .4,3(13).  .4  J2l  (62  -  .42.31(3).  U3). 
In  [14].  numerical  experiments  were  performed  with  these  and  other  preconditioners. 


6.  A  new  class  of  banded,  row-sum  preserving  preconditioners 

We  now  present  a  new  family  of  preconditioned  for  the  capacitance  matrix  C.  These  pre- 
conditioners  are  motivated  by  the  empirical  observation  that  the  elements  of  the  matrix  C  decay 
away  from  the  main  diagonal.  It  is.  therefore,  reasonable  to  consider  ^--diagonal  approximations 
to  C.  It  would  not.  however,  be  efficient  to  compute  the  elements  of  C  in  order  to  do  this.  We 
now  present  a  method  for  computing  a  ^--diagonal  approximation  to  C  without  requiring  the  com¬ 
putation  of  C  explicitly.  The  idea  is  motivated  by  sparse  Jacobian  evaluation  techniques  [9].  For 
example,  for  the  case  k  =  3.  the  approximant  M  to  C  can  be  computed  in  compact  form  by  eval¬ 
uating  the  three  products  Ck,.  /  =  1.  2.3,  where  u-%  =  (1,0,0. 1.0. . . . )T ,  *23  —  (0.  1,0.  0. 1. . .  ,)T  and 
us  =  (0.0.  1.0.0. . .  d7".  The  motivation  is  clear,  for  if  C  were  indeed  tridiagonal,  (k  =  3).  then  all  of 
its  nonzero  elements  can  be  found  in  the  three  vectors  Cu,.  1  =  1. 2.3.  Note  that  the  computation  of 
each  product  Cu,  involves  solving  one  problem  on  each  subdomain  with  u,  as  boundary  condition 
on  the  interface. 

The  generalization  to  other  values  of  k  is  obvious.  Moreover,  it  can  be  easily  verified  that  the 
matrix  M  computed  this  way  preserves  the  row-sums  of  ('.  The  case  k  —  1.  however,  deserves 
special  mention.  The  method  described  above  would  compute  a  diagonal  approximation  to  C .  with 


r  <r  [  -r  -  .  w" 


~  irv*r%  *r ,  ,  *  ,  -  <  —  *»  -  ■_ 


diagonal  entries  given  by  Ce.  where  e  -  (1,1 . I)7".  However,  since  the  first  term  A33  in  the 

definition  of  C  in  (2.4)  is  already  known  explicitly  (and  it  is  tridiagonal),  it  is  only  necessary  to 
apply  the  above  approximation  procedure  to  the  last  two  terms  in  (2.4).  The  resulting  matrix  A/  is 
thus  tridiagonal,  namely.  A33  with  the  diagonal  entries  modified  in  such  a  way  that  the  row  sums 
of  C  are  preserved.  Viewed  this  way,  the  case  k  —  1  is  similar  in  spirit  to  the  Dupont-Ivendall- 
Rachford  procedure  [  1 1]  for  obtaining  an  easily  invertible  banded  approximant  for  C .  This  special 
procedure  for  the  case  k  =  1  was  sugested  independently  by  Eisenstat  [12].  See  [  14]  for  numerical 
experiments  with  this  class  of  preconditioners. 

In  general,  for  a  ^-diagonal  approximation  to  C.  k  problems  on  each  subdomain  must  be 
solved,  which  may  seem  prohibitively  expensive  except  for  small  values  of  k.  However,  the  main 
advantage  of  this  family  of  preconditioners  is  that  they  are  less  dependent  on  special  properties 
(e.g.  eigeustructures)  of  the  differential  operator  underlying  A.  Moreover,  for  nonlinear  problems 
where  a  Newton  type  outer  iteration  may  be  involved,  one  preconditioner  can  be  reused  several 
times  and  the  cost  of  computing  it  can  be  amortized  over  the  overall  iteration.  Further  details  will 
be  reported  elsewhere. 
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